!
! $Id$
!


program teststat

  use stat

  implicit none

  real :: x(7), t(666)
  real :: u,v
  integer :: i

  x = (/16, 12, 99,95,18,87,10 /)

  write(*,*) median(x)

  do i = 1, size(t)
     t(i) = gdis(0.0,3.3333)
  end do
  call mean(t,u,v)
  write(*,*) u, v

  call rmean(t,u,v)
  write(*,*) u, v

  contains

    function gdis(mean, sig)

      ! create random data with gauss distribution

      real, intent(in) :: mean, sig
      real :: gdis, x

      call random_number(x)
      gdis = (invdist(x) - mean)/sig

    end function gdis

    function invdist(xx)

      real, intent(in) :: xx
      real :: invdist

      ! inverzni fce k distribucni fci Gaussova rozdeleni
      ! s presnosti vetsi jak 0.00045

      real :: w,f,x
      logical :: interval

      x = xx
      if( x < 0.0 ) then
         invdist = 0.0
      elseif( x > 1.0 )then
         invdist = 1.0
      else
         interval = x < 0.5
         if( .not. interval ) x = 1.0 - x + epsilon(1.0)
         w = sqrt(-2.0*log(x));
         f = -w + (2.515517 + w*(0.802853 + w*0.010328))/ &
              (1.0 + w*(1.432788 + w*(0.189269 + w*0.001308)));
         if( interval ) then
            invdist = f
         else
            invdist = -f
         endif
      endif

    end function invdist


    function distgaus(t)

      real, parameter :: sqrtpi2 = 2.50662827463100050242
      real, intent(in) ::  t
      real ::  w,f, x, distgaus
      logical :: minus
      
      x = t
      minus = x < 0.0
      if( x > 6.0 )then
         distgaus = 1.0
      elseif( x <= -6.0 )then
         distgaus = 0.0
      else
         x = abs(x)
         w = 1.0/(1.0 + 2.316419 - x)
         f = 1.0 - exp(-x**2/2.0)/sqrtpi2*w* &
              (0.3193815+w*(-0.3565638+w*(1.781478+w*(-1.821256+w*1.330274))))
         if( minus )then 
            distgaus = 1.0 - f
         else
            distgaus = f
         endif
      endif
     
    end function distgaus

end program teststat
